APIM - Control - Pushing Separate

knitr::opts_chunk$set(echo = TRUE)

cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)

system("shutdown /a")

## [1] 1116

# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
report_hurdle = TRUE
do_priorsense = FALSE

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 2000 # 10'000 per chain to achieve 40'000
warmup = 1000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = as.character(iterations)

df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = FALSE, use_mi = use_mi)[[1]]

Constructing scales reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.3468 0.0000 5.0000 2323

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'planPlan',
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')

# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Actionplan",
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')

rows_to_pack <- list(
  "Within-Person Effects" = c(2,10),
  "Between-Person Effects" = c(11,17),
  "Random Effects" = c(18, 24), 
  "Additional Parameters" = c(25,25)
  )


rows_to_pack_ordinal <- list(
  "Intercepts" = c(1,6),
  "Within-Person Effects" = c(2+5,10+5),
  "Between-Person Effects" = c(11+5,17+5),
  "Random Effects" = c(18+5, 24+5), 
  "Additional Parameters" = c(25+5,25+6)
  )

HURDLE MODELS

# For indistinguishable Dyads
model_rows_fixed_hu <- c(
    'Intercept', 
    'hu_Intercept',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'planPlan',
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb',
    
    # HURDLE MODEL
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'hu_persuasion_self_cw', 
    'hu_persuasion_partner_cw', 
    'hu_pressure_self_cw', 
    'hu_pressure_partner_cw', 
    'hu_pushing_self_cw', 
    'hu_pushing_partner_cw', 
    'hu_day', 
    'hu_planPlan',
    'hu_weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'hu_persuasion_self_cb',
    'hu_persuasion_partner_cb',
    'hu_pressure_self_cb',
    'hu_pressure_partner_cb',
    'hu_pushing_self_cb',
    'hu_pushing_partner_cb',
    'hu_weartime_self_cb'
  )

model_rows_fixed_hu_ordinal <- c(
  model_rows_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed_hu[3:length(model_rows_fixed_hu)]
)


model_rows_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(hu_Intercept)',
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # HURDLE
  'sd(hu_persuasion_self_cw)',
  'sd(hu_persuasion_partner_cw)',
  'sd(hu_pressure_self_cw)',
  'sd(hu_pressure_partner_cw)',
  'sd(hu_pushing_self_cw)',
  'sd(hu_pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_hu_ordinal <- c(model_rows_random_hu,'disc')

# For indistinguishable Dyads
model_rownames_fixed_hu <- c(
    "Intercept", 
    "Hurdle Intercept",
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Actionplan",
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime",
    
    # HURDLE
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Hu Daily persuasion experienced", 
    "Hu Daily persuasion utilized (partner's view)", # OR partner received
    "Hu Daily pressure experienced", 
    "Hu Daily pressure utilized (partner's view)", 
    "Hu Daily pushing experienced", 
    "Hu Daily pushing utilized (partner's view)", 
    "Hu Day", 
    "Hu Actionplan",
    "Hu Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Hu Mean persuasion experienced", 
    "Hu Mean persuasion utilized (partner's view)", 
    "Hu Mean pressure experienced", 
    "Hu Mean pressure utilized (partner's view)", 
    "Hu Mean pushing experienced", 
    "Hu Mean pushing utilized (partner's view)", 
    "Hu Mean weartime"
  )



model_rownames_fixed_hu_ordinal <- c(
  model_rownames_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed_hu[3:length(model_rownames_fixed_hu)]
)



model_rownames_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Hurdle Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  
  # Hurdle
  "sd(Hu Daily persuasion experienced)", 
  "sd(Hu Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Hu Daily pressure experienced)", 
  "sd(Hu Daily pressure utilized (partner's view))", 
  "sd(Hu Daily pushing experienced)", 
  "sd(Hu Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_hu_ordinal <- c(model_rownames_random_hu,'disc')

rows_to_pack_hu <- list(
  "Conditional Within-Person Effects" = c(3,11),
  "Conditional Between-Person Effects" = c(12,18),
  
  "Hurdle Within-Person Effects" = c(19,27),
  "Hurdle Between-Person Effects" = c(28,34),
  
  "Random Effects" = c(35, 48), 
  "Additional Parameters" = c(49,49)
  )

rows_to_pack_hu_ordinal <- list(
  "Intercepts" = c(1,7),
  "Conditional Within-Person Effects" = c(3+5,11+5),
  "Conditional Between-Person Effects" = c(12+5,18+5),
  
  "Hurdle Within-Person Effects" = c(19+5,27+5),
  "Hurdle Between-Person Effects" = c(28+5,34+5),
  
  "Random Effects" = c(35+5, 48+5), 
  "Additional Parameters" = c(49+5,49+6)
  )

Self-Reported MVPA

range(df_double$pa_sub, na.rm = T) 

## [1]   0 720

hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

WITHOUT PUSHING

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    
    plan +
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    
    plan +
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_Nopushing", suffix))
  #, file_refit = 'always'
)

## Warning: Rows containing NAs were excluded from the model.

pa_sub_digest <- digest::digest(pa_sub)

# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
  pa_sub, 
  stats_to_report = c('pd','ROPE'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 

## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.

# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub_continuous <- summarize_brms(
  pa_sub, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_continuous,
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 

## Warning in summarize_brms(pa_sub, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.

# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous

columns_to_replace <- c("ROPE", "inside ROPE")

summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <- 
  summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]

# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack_hu)

	exp(Est.)	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	39.09***	2.57	[34.11, 44.31]	1.000	[0.92, 1.08]	0.000	1.000	1519	2186
Hurdle Intercept	0.27***	0.04	[ 0.20, 0.36]	1.000	[0.84, 1.20]	0.000	1.002	1267	1972
Conditional Within-Person Effects
Daily persuasion experienced	1.02	0.03	[ 0.97, 1.08]	0.828	[0.92, 1.08]	0.980	1.000	2047	2905
Daily persuasion utilized (partner’s view)	1.02	0.02	[ 0.98, 1.07]	0.866	[0.92, 1.08]	0.992	1.000	2764	3010
Daily pressure experienced	0.91*	0.04	[ 0.82, 1.00]	0.978	[0.92, 1.08]	0.329	1.001	4179	2899
Daily pressure utilized (partner’s view)	0.93	0.04	[ 0.85, 1.02]	0.944	[0.92, 1.08]	0.590	1.000	4011	2995
Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Day	0.96	0.06	[ 0.86, 1.09]	0.723	[0.92, 1.08]	0.722	1.001	5576	2644
Actionplan	1.38***	0.06	[ 1.27, 1.50]	1.000	[0.92, 1.08]	0.000	1.002	6770	3219
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Conditional Between-Person Effects
Mean persuasion experienced	1.10	0.15	[ 0.83, 1.44]	0.757	[0.92, 1.08]	0.335	1.001	924	1425
Mean persuasion utilized (partner’s view)	1.07	0.15	[ 0.82, 1.41]	0.692	[0.92, 1.08]	0.382	1.003	941	1795
Mean pressure experienced	1.15	0.19	[ 0.84, 1.60]	0.801	[0.92, 1.08]	0.268	1.001	1275	2139
Mean pressure utilized (partner’s view)	0.96	0.16	[ 0.69, 1.33]	0.599	[0.92, 1.08]	0.347	1.002	1199	1938
Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced	1.55***	0.11	[ 1.36, 1.79]	1.000	[0.84, 1.20]	0.000	1.001	2972	2563
Hu Daily persuasion utilized (partner’s view)	1.37***	0.08	[ 1.23, 1.57]	1.000	[0.84, 1.20]	0.007	1.001	4448	2648
Hu Daily pressure experienced	0.99	0.13	[ 0.74, 1.32]	0.545	[0.84, 1.20]	0.796	1.001	4751	2981
Hu Daily pressure utilized (partner’s view)	1.59**	0.28	[ 1.14, 2.40]	0.997	[0.84, 1.20]	0.048	1.001	4149	2710
Hu Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Day	0.85	0.12	[ 0.64, 1.12]	0.868	[0.84, 1.20]	0.538	1.004	7877	3012
Hu Actionplan	9.74***	0.88	[ 8.17, 11.74]	1.000	[0.84, 1.20]	0.000	1.005	8410	3039
Hu Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced	1.53	0.48	[ 0.80, 2.84]	0.902	[0.84, 1.20]	0.194	1.002	1143	1885
Hu Mean persuasion utilized (partner’s view)	1.58	0.50	[ 0.85, 2.94]	0.923	[0.84, 1.20]	0.169	1.003	1108	1848
Hu Mean pressure experienced	0.34**	0.13	[ 0.16, 0.73]	0.997	[0.84, 1.20]	0.010	1.002	1691	2335
Hu Mean pressure utilized (partner’s view)	0.56	0.22	[ 0.26, 1.20]	0.930	[0.84, 1.20]	0.132	1.002	1612	2154
Hu Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.30	0.04	[0.23, 0.40]	NA	NA	NA	1.003	1023	2027
sd(Hurdle Intercept)	0.73	0.10	[0.55, 0.97]	NA	NA	NA	1.001	1448	2056
sd(Daily persuasion experienced)	0.12	0.02	[0.08, 0.17]	NA	NA	NA	1.001	2526	2989
sd(Daily persuasion utilized (partner’s view))	0.08	0.02	[0.05, 0.13]	NA	NA	NA	1.000	2678	2238
sd(Daily pressure experienced)	0.07	0.06	[0.00, 0.24]	NA	NA	NA	1.005	1479	1735
sd(Daily pressure utilized (partner’s view))	0.07	0.05	[0.00, 0.19]	NA	NA	NA	1.002	2054	2288
sd(Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily persuasion experienced)	0.23	0.08	[0.07, 0.41]	NA	NA	NA	1.002	1302	1082
sd(Hu Daily persuasion utilized (partner’s view))	0.16	0.09	[0.01, 0.33]	NA	NA	NA	1.002	1146	997
sd(Hu Daily pressure experienced)	0.17	0.15	[0.01, 0.66]	NA	NA	NA	1.003	1846	2077
sd(Hu Daily pressure utilized (partner’s view))	0.22	0.20	[0.01, 0.86]	NA	NA	NA	1.001	1680	2166
sd(Hu Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
Additional Parameters
sigma	0.68	0.01	[0.66, 0.70]	NA	NA	NA	1.001	7369	2900

# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()

## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.77), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.79), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.78), b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.73), b_hu_pressure_self_cb and
##   b_hu_persuasion_self_cb (r = 0.76), b_hu_pressure_partner_cb and
##   b_hu_persuasion_partner_cb (r = 0.76). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

Device Based MVPA

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    
    plan +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_nopushing", suffix))
)

## Warning: Rows containing NAs were excluded from the model.

pa_obj_log_digest <- digest::digest(pa_obj_log)

# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summarize_brms(
  pa_obj_log, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)

## Warning in summarize_brms(pa_obj_log, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.

	exp(Est.)	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	111.65***	5.82	[101.28, 124.78]	1.000	[0.94, 1.07]	0.000	1.007	717	1643
Within-Person Effects
Daily persuasion experienced	1.03	0.02	[ 1.00, 1.06]	0.973	[0.94, 1.07]	0.990	1.000	2672	3035
Daily persuasion utilized (partner’s view)	1.02	0.02	[ 0.99, 1.05]	0.892	[0.94, 1.07]	0.997	1.000	3291	3157
Daily pressure experienced	0.95	0.03	[ 0.89, 1.02]	0.926	[0.94, 1.07]	0.685	1.001	5011	2874
Daily pressure utilized (partner’s view)	0.98	0.03	[ 0.92, 1.05]	0.691	[0.94, 1.07]	0.930	1.001	4611	3115
Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Day	0.97	0.04	[ 0.90, 1.04]	0.777	[0.94, 1.07]	0.837	1.003	7630	2798
Actionplan	1.09***	0.02	[ 1.05, 1.14]	1.000	[0.94, 1.07]	0.144	1.000	8010	2892
Daily weartime	1.00***	0.00	[ 1.00, 1.00]	1.000	[0.94, 1.07]	1.000	1.000	4046	2780
Between-Person Effects
Mean persuasion experienced	1.10	0.14	[ 0.86, 1.42]	0.787	[0.94, 1.07]	0.306	1.005	796	1443
Mean persuasion utilized (partner’s view)	1.02	0.13	[ 0.80, 1.32]	0.573	[0.94, 1.07]	0.392	1.006	784	1435
Mean pressure experienced	0.91	0.13	[ 0.69, 1.18]	0.749	[0.94, 1.07]	0.294	1.004	982	1686
Mean pressure utilized (partner’s view)	1.04	0.14	[ 0.79, 1.35]	0.595	[0.94, 1.07]	0.351	1.008	842	1595
Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean weartime	1.00	0.00	[ 1.00, 1.00]	0.944	[0.94, 1.07]	1.000	1.001	4896	3276
Random Effects
sd(Intercept)	0.29	0.04	[0.23, 0.38]	NA	NA	NA	1.001	876	1677
sd(Daily persuasion experienced)	0.05	0.01	[0.02, 0.08]	NA	NA	NA	1.002	2216	2269
sd(Daily persuasion utilized (partner’s view))	0.05	0.02	[0.02, 0.09]	NA	NA	NA	1.001	1646	1915
sd(Daily pressure experienced)	0.04	0.04	[0.00, 0.15]	NA	NA	NA	1.000	1899	2148
sd(Daily pressure utilized (partner’s view))	0.04	0.03	[0.00, 0.11]	NA	NA	NA	1.000	2228	2102
sd(Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
Additional Parameters
sigma	0.58	0.01	[0.56, 0.59]	NA	NA	NA	1.002	5669	3240

plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()

## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.86), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.84), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.78), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.78), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.86). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    plan +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss_nopushing", suffix))
)

## Warning: Rows containing NAs were excluded from the model.

mood_gauss_digest <- digest::digest(mood_gauss)

summarize_brms(
  mood_gauss, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)

	Est.	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	3.66***	0.11	[ 3.45, 3.86]	1.000	[-0.11, 0.11]	0.000	1.002	349	743
Within-Person Effects
Daily persuasion experienced	0.00	0.02	[-0.04, 0.05]	0.565	[-0.11, 0.11]	1.000	1.001	3204	2834
Daily persuasion utilized (partner’s view)	0.03	0.02	[-0.01, 0.08]	0.912	[-0.11, 0.11]	0.999	1.000	2517	2838
Daily pressure experienced	-0.03	0.05	[-0.13, 0.07]	0.721	[-0.11, 0.11]	0.944	1.000	3693	3168
Daily pressure utilized (partner’s view)	0.00	0.05	[-0.10, 0.11]	0.535	[-0.11, 0.11]	0.967	1.001	3267	2612
Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Day	0.25***	0.06	[ 0.14, 0.37]	1.000	[-0.11, 0.11]	0.009	1.000	6114	2852
Actionplan	0.10**	0.03	[ 0.03, 0.17]	0.999	[-0.11, 0.11]	0.617	1.002	6449	3035
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean persuasion experienced	0.40	0.25	[-0.10, 0.88]	0.946	[-0.11, 0.11]	0.108	1.003	512	987
Mean persuasion utilized (partner’s view)	0.31	0.25	[-0.18, 0.79]	0.899	[-0.11, 0.11]	0.178	1.003	518	889
Mean pressure experienced	-0.32	0.27	[-0.85, 0.20]	0.890	[-0.11, 0.11]	0.169	1.004	575	1114
Mean pressure utilized (partner’s view)	-0.26	0.27	[-0.78, 0.27]	0.830	[-0.11, 0.11]	0.219	1.003	581	1179
Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.59	0.08	[0.47, 0.78]	NA	NA	NA	1.008	658	1266
sd(Daily persuasion experienced)	0.05	0.03	[0.00, 0.10]	NA	NA	NA	1.003	906	1079
sd(Daily persuasion utilized (partner’s view))	0.08	0.03	[0.03, 0.13]	NA	NA	NA	1.001	1496	1601
sd(Daily pressure experienced)	0.06	0.06	[0.00, 0.22]	NA	NA	NA	1.000	1716	1674
sd(Daily pressure utilized (partner’s view))	0.07	0.06	[0.00, 0.24]	NA	NA	NA	1.002	1396	1589
sd(Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
Additional Parameters
sigma	0.96	0.01	[0.94, 0.98]	NA	NA	NA	1.000	6259	3060

plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()

## Possible multicollinearity between b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.9), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.84), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.85), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.89), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.8). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

Reactance

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    plan +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal_nopushing", suffix))
)

## Warning: Rows containing NAs were excluded from the model.

reactance_ordinal_digest <- digest::digest(reactance_ordinal)

summarize_brms(
  reactance_ordinal, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)

	OR	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercepts
Intercept	NA	NA	NA	NA	NA	NA	NA	NA	NA
Intercept[1]	3.58***	1.04	[ 2.03, 6.38]	1.000	[0.84, 1.20]	0.000	1.000	3499	2918
Intercept[2]	7.49***	2.25	[ 4.22, 13.84]	1.000	[0.84, 1.20]	0.000	0.999	3448	2962
Intercept[3]	20.32***	6.43	[ 11.09, 39.16]	1.000	[0.84, 1.20]	0.000	1.000	3542	3144
Intercept[4]	87.96***	32.35	[ 42.93, 187.99]	1.000	[0.84, 1.20]	0.000	1.001	4026	3073
Intercept[5]	2856.04***	1828.76	[894.50, 11588.88]	1.000	[0.84, 1.20]	0.000	1.000	5154	3183
Within-Person Effects
Daily persuasion experienced	0.85*	0.07	[ 0.72, 0.99]	0.983	[0.84, 1.20]	0.593	1.001	3783	3010
Daily persuasion utilized (partner’s view)	1.01	0.09	[ 0.84, 1.20]	0.548	[0.84, 1.20]	0.949	1.001	3102	2913
Daily pressure experienced	2.01**	0.35	[ 1.33, 2.82]	0.998	[0.84, 1.20]	0.009	1.001	2513	2764
Daily pressure utilized (partner’s view)	1.15	0.23	[ 0.73, 1.84]	0.761	[0.84, 1.20]	0.516	1.001	2456	2104
Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Day	1.64	0.56	[ 0.84, 3.18]	0.929	[0.84, 1.20]	0.154	1.000	5945	3187
Actionplan	1.03	0.21	[ 0.70, 1.50]	0.551	[0.84, 1.20]	0.618	1.000	4672	2843
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean persuasion experienced	1.08	0.51	[ 0.42, 2.80]	0.572	[0.84, 1.20]	0.302	1.004	1403	2385
Mean persuasion utilized (partner’s view)	0.78	0.41	[ 0.30, 2.20]	0.677	[0.84, 1.20]	0.238	1.005	1606	2410
Mean pressure experienced	3.32*	1.70	[ 1.19, 9.62]	0.985	[0.84, 1.20]	0.019	1.004	1692	2235
Mean pressure utilized (partner’s view)	1.23	0.68	[ 0.40, 3.42]	0.647	[0.84, 1.20]	0.236	1.004	1685	2580
Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.78	0.19	[0.47, 1.22]	NA	NA	NA	1.002	1599	2400
sd(Daily persuasion experienced)	0.17	0.13	[0.01, 0.43]	NA	NA	NA	1.006	638	1049
sd(Daily persuasion utilized (partner’s view))	0.20	0.12	[0.02, 0.48]	NA	NA	NA	1.004	1126	1063
sd(Daily pressure experienced)	0.47	0.22	[0.07, 1.01]	NA	NA	NA	1.003	1001	847
sd(Daily pressure utilized (partner’s view))	0.29	0.29	[0.01, 1.41]	NA	NA	NA	1.001	1145	2311
sd(Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
Additional Parameters
sigma	NA	NA	NA	NA	NA	NA	NA	NA	NA
disc	1.00	0.00	[1.00, 1.00]	NA	NA	NA	NA	NA	NA

plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()

## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.8), b_Intercept[4] and b_Intercept[3] (r = 0.86), b_pressure_self_cb
##   and b_persuasion_self_cb (r = 0.8), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.84). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conditional_spaghetti(
  reactance_ordinal, 
  effects = c('persuasion_self_cw', 'pressure_self_cw')
  , group_var = 'coupleID'
  #, n_groups = 15
  , plot_full_range = T
)

\(persuasion_self_cw <img src="03_SensitivityPushingSeparate_files/figure-html/report_reactance_ordinal-3.png" width="2400" />\)pressure_self_cw

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    plan +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance_nopushing_", suffix))
  #, file_refit = 'always'
)

## Warning: Rows containing NAs were excluded from the model.

is_reactance_digest <- digest::digest(is_reactance)

summarize_brms(
  is_reactance, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)

	OR	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	0.31***	0.10	[0.16, 0.57]	1.000	[0.83, 1.20]	0.000	1.001	4202	3377
Within-Person Effects
Daily persuasion experienced	0.85	0.08	[0.71, 1.02]	0.957	[0.83, 1.20]	0.605	1.000	3643	3198
Daily persuasion utilized (partner’s view)	1.10	0.14	[0.86, 1.46]	0.793	[0.83, 1.20]	0.719	1.000	3319	3084
Daily pressure experienced	2.09*	0.53	[1.25, 4.16]	0.994	[0.83, 1.20]	0.020	1.001	2178	1389
Daily pressure utilized (partner’s view)	1.33	0.50	[0.63, 4.00]	0.786	[0.83, 1.20]	0.287	1.003	1911	1442
Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Day	1.84	0.68	[0.90, 3.72]	0.957	[0.83, 1.20]	0.107	1.001	6332	2915
Actionplan	1.07	0.25	[0.67, 1.74]	0.622	[0.83, 1.20]	0.542	0.999	5977	3011
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean persuasion experienced	1.91	1.12	[0.63, 6.06]	0.871	[0.83, 1.20]	0.142	1.002	1488	1842
Mean persuasion utilized (partner’s view)	1.07	0.66	[0.33, 3.87]	0.539	[0.83, 1.20]	0.226	1.001	1732	2699
Mean pressure experienced	17.22***	15.48	[3.14, 113.96]	1.000	[0.83, 1.20]	0.001	1.001	2115	2134
Mean pressure utilized (partner’s view)	1.15	1.09	[0.17, 7.83]	0.556	[0.83, 1.20]	0.146	1.002	2328	3018
Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	1.21	0.24	[0.81, 1.76]	NA	NA	NA	1.004	1669	2622
sd(Daily persuasion experienced)	0.19	0.13	[0.01, 0.50]	NA	NA	NA	1.008	908	1526
sd(Daily persuasion utilized (partner’s view))	0.43	0.17	[0.12, 0.84]	NA	NA	NA	1.002	1349	988
sd(Daily pressure experienced)	0.77	0.53	[0.05, 2.08]	NA	NA	NA	1.004	775	1327
sd(Daily pressure utilized (partner’s view))	0.74	0.66	[0.03, 2.59]	NA	NA	NA	1.006	1088	1559
sd(Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
Additional Parameters
sigma	NA	NA	NA	NA	NA	NA	NA	NA	NA

plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

Report All Models

process_model_component <- function(obj) {
  # Convert to string, modify, and evaluate
  name <- deparse(substitute(obj))
  if (report_hurdle) name <- paste0(name, '_hu')
  if (report_ordinal) name <- paste0(name, '_ordinal')
  return(get(name, envir = parent.frame()))
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  reactance_ordinal,
  is_reactance,
  
  stats_to_report = c('CI', 'pd'),
  
  model_rows_random = process_model_component(model_rows_random),
  model_rows_fixed = process_model_component(model_rows_fixed),
  model_rownames_random = process_model_component(model_rownames_random),
  model_rownames_fixed = process_model_component(model_rownames_fixed)
) 

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”

# pretty printing

summary_all_models <- all_models %>%
  print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),  
      "Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5), 
      "Mood Gaussian" = (length(all_models) / 5),
      "Reactance Ordinal" = (length(all_models) / 5),
      "Reactance Dichotome" = (length(all_models) / 5))
  )


export_xlsx(
  summary_all_models, 
  rows_to_pack = process_model_component(rows_to_pack),
  file.path("Output", "AllModelsNoPushing.xlsx"), 
  merge_option = 'header', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

## 
## Attaching package: 'rvest'

## The following object is masked from 'package:readr':
## 
##     guess_encoding

print(summary_all_models)

	Subjective MVPA Hurdle Lognormal			Device-Based MVPA Log (Gaussian)			Mood Gaussian			Reactance Ordinal			Reactance Dichotome
	exp(Est.) pa_sub	95% CI pa_sub	pd pa_sub	exp(Est.) pa_obj_log	95% CI pa_obj_log	pd pa_obj_log	Est. mood_gauss	95% CI mood_gauss	pd mood_gauss	OR reactance_ordinal	95% CI reactance_ordinal	pd reactance_ordinal	OR is_reactance	95% CI is_reactance	pd is_reactance
Intercept	39.09***	[34.11, 44.31]	1.000	111.65***	[101.28, 124.78]	1.000	3.66***	[ 3.45, 3.86]	1.000	NA	NA	NA	0.31***	[0.16, 0.57]	1.000
Hurdle Intercept	0.27***	[ 0.20, 0.36]	1.000	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Conditional Within-Person Effects
Daily persuasion experienced	1.02	[ 0.97, 1.08]	0.828	1.03	[ 1.00, 1.06]	0.973	0.00	[-0.04, 0.05]	0.565	0.85*	[ 0.72, 0.99]	0.983	0.85	[0.71, 1.02]	0.957
Daily persuasion utilized (partner’s view)	1.02	[ 0.98, 1.07]	0.866	1.02	[ 0.99, 1.05]	0.892	0.03	[-0.01, 0.08]	0.912	1.01	[ 0.84, 1.20]	0.548	1.10	[0.86, 1.46]	0.793
Daily pressure experienced	0.91*	[ 0.82, 1.00]	0.978	0.95	[ 0.89, 1.02]	0.926	-0.03	[-0.13, 0.07]	0.721	2.01**	[ 1.33, 2.82]	0.998	2.09*	[1.25, 4.16]	0.994
Daily pressure utilized (partner’s view)	0.93	[ 0.85, 1.02]	0.944	0.98	[ 0.92, 1.05]	0.691	0.00	[-0.10, 0.11]	0.535	1.15	[ 0.73, 1.84]	0.761	1.33	[0.63, 4.00]	0.786
Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Day	0.96	[ 0.86, 1.09]	0.723	0.97	[ 0.90, 1.04]	0.777	0.25***	[ 0.14, 0.37]	1.000	1.64	[ 0.84, 3.18]	0.929	1.84	[0.90, 3.72]	0.957
Actionplan	1.38***	[ 1.27, 1.50]	1.000	1.09***	[ 1.05, 1.14]	1.000	0.10**	[ 0.03, 0.17]	0.999	1.03	[ 0.70, 1.50]	0.551	1.07	[0.67, 1.74]	0.622
Daily weartime	NA	NA	NA	1.00***	[ 1.00, 1.00]	1.000	NA	NA	NA	NA	NA	NA	NA	NA	NA
Conditional Between-Person Effects
Mean persuasion experienced	1.10	[ 0.83, 1.44]	0.757	1.10	[ 0.86, 1.42]	0.787	0.40	[-0.10, 0.88]	0.946	1.08	[ 0.42, 2.80]	0.572	1.91	[0.63, 6.06]	0.871
Mean persuasion utilized (partner’s view)	1.07	[ 0.82, 1.41]	0.692	1.02	[ 0.80, 1.32]	0.573	0.31	[-0.18, 0.79]	0.899	0.78	[ 0.30, 2.20]	0.677	1.07	[0.33, 3.87]	0.539
Mean pressure experienced	1.15	[ 0.84, 1.60]	0.801	0.91	[ 0.69, 1.18]	0.749	-0.32	[-0.85, 0.20]	0.890	3.32*	[ 1.19, 9.62]	0.985	17.22***	[3.14, 113.96]	1.000
Mean pressure utilized (partner’s view)	0.96	[ 0.69, 1.33]	0.599	1.04	[ 0.79, 1.35]	0.595	-0.26	[-0.78, 0.27]	0.830	1.23	[ 0.40, 3.42]	0.647	1.15	[0.17, 7.83]	0.556
Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean weartime	NA	NA	NA	1.00	[ 1.00, 1.00]	0.944	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced	1.55***	[ 1.36, 1.79]	1.000	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily persuasion utilized (partner’s view)	1.37***	[ 1.23, 1.57]	1.000	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pressure experienced	0.99	[ 0.74, 1.32]	0.545	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pressure utilized (partner’s view)	1.59**	[ 1.14, 2.40]	0.997	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Day	0.85	[ 0.64, 1.12]	0.868	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Actionplan	9.74***	[ 8.17, 11.74]	1.000	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced	1.53	[ 0.80, 2.84]	0.902	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean persuasion utilized (partner’s view)	1.58	[ 0.85, 2.94]	0.923	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pressure experienced	0.34**	[ 0.16, 0.73]	0.997	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pressure utilized (partner’s view)	0.56	[ 0.26, 1.20]	0.930	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pushing experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pushing utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.30	[0.23, 0.40]	NA	0.29	[0.23, 0.38]	NA	0.59	[0.47, 0.78]	NA	0.78	[0.47, 1.22]	NA	1.21	[0.81, 1.76]	NA
sd(Hurdle Intercept)	0.73	[0.55, 0.97]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion experienced)	0.12	[0.08, 0.17]	NA	0.05	[0.02, 0.08]	NA	0.05	[0.00, 0.10]	NA	0.17	[0.01, 0.43]	NA	0.19	[0.01, 0.50]	NA
sd(Daily persuasion utilized (partner’s view))	0.08	[0.05, 0.13]	NA	0.05	[0.02, 0.09]	NA	0.08	[0.03, 0.13]	NA	0.20	[0.02, 0.48]	NA	0.43	[0.12, 0.84]	NA
sd(Daily pressure experienced)	0.07	[0.00, 0.24]	NA	0.04	[0.00, 0.15]	NA	0.06	[0.00, 0.22]	NA	0.47	[0.07, 1.01]	NA	0.77	[0.05, 2.08]	NA
sd(Daily pressure utilized (partner’s view))	0.07	[0.00, 0.19]	NA	0.04	[0.00, 0.11]	NA	0.07	[0.00, 0.24]	NA	0.29	[0.01, 1.41]	NA	0.74	[0.03, 2.59]	NA
sd(Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily persuasion experienced)	0.23	[0.07, 0.41]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily persuasion utilized (partner’s view))	0.16	[0.01, 0.33]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pressure experienced)	0.17	[0.01, 0.66]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pressure utilized (partner’s view))	0.22	[0.01, 0.86]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pushing experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pushing utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Additional Parameters
sigma	0.68	[0.66, 0.70]	NA	0.58	[0.56, 0.59]	NA	0.96	[0.94, 0.98]	NA	NA	NA	NA	NA	NA	NA

ONLY PUSHING

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan +
    day + 
    
    # Random effects
    (1 + pushing_self_cw + pushing_partner_cw | coupleID),
  
  hu = ~ pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan +
    day + 
    
    # Random effects
    (1 + pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_onlypushing", suffix))
  #, file_refit = 'always'
)

## Warning: Rows containing NAs were excluded from the model.

pa_sub_digest <- digest::digest(pa_sub)

# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
  pa_sub, 
  stats_to_report = c('pd','ROPE'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 

## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.

# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub_continuous <- summarize_brms(
  pa_sub, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_continuous,
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 

## Warning in summarize_brms(pa_sub, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.

# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous

columns_to_replace <- c("ROPE", "inside ROPE")

summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <- 
  summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]

# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack_hu)

	exp(Est.)	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	57.28***	4.53	[49.10, 66.97]	1.000	[0.93, 1.08]	0.000	1.002	1327	2040
Hurdle Intercept	3.79***	0.98	[ 2.27, 6.26]	1.000	[0.84, 1.20]	0.000	1.005	1088	2463
Conditional Within-Person Effects
Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing experienced	0.97	0.03	[ 0.91, 1.03]	0.856	[0.93, 1.08]	0.931	1.002	4427	2662
Daily pushing utilized (partner’s view)	0.97	0.03	[ 0.90, 1.03]	0.847	[0.93, 1.08]	0.907	1.001	3679	2873
Day	0.91	0.07	[ 0.79, 1.06]	0.894	[0.93, 1.08]	0.398	1.000	8744	2600
Actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Conditional Between-Person Effects
Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing experienced	1.32*	0.18	[ 1.00, 1.71]	0.977	[0.93, 1.08]	0.060	1.005	1217	2176
Mean pushing utilized (partner’s view)	1.18	0.15	[ 0.90, 1.53]	0.883	[0.93, 1.08]	0.222	1.005	1265	2076
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pushing experienced	1.55**	0.23	[ 1.17, 2.24]	0.998	[0.84, 1.20]	0.035	1.002	3127	2531
Hu Daily pushing utilized (partner’s view)	1.84***	0.30	[ 1.38, 2.68]	1.000	[0.84, 1.20]	0.002	1.001	2766	2989
Hu Day	0.61	0.15	[ 0.37, 1.00]	0.974	[0.84, 1.20]	0.100	1.003	7632	2768
Hu Actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pushing experienced	0.46**	0.14	[ 0.26, 0.81]	0.997	[0.84, 1.20]	0.018	1.001	2675	2824
Hu Mean pushing utilized (partner’s view)	0.62	0.18	[ 0.33, 1.08]	0.952	[0.84, 1.20]	0.136	1.001	2568	2699
Hu Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.37	0.06	[0.29, 0.50]	NA	NA	NA	1.002	1283	2370
sd(Hurdle Intercept)	1.25	0.19	[0.93, 1.74]	NA	NA	NA	1.002	1448	2426
sd(Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing experienced)	0.06	0.04	[0.00, 0.16]	NA	NA	NA	1.001	870	1486
sd(Daily pushing utilized (partner’s view))	0.09	0.04	[0.01, 0.18]	NA	NA	NA	1.004	977	1142
sd(Hu Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pushing experienced)	0.46	0.24	[0.07, 1.02]	NA	NA	NA	1.000	947	1689
sd(Hu Daily pushing utilized (partner’s view))	0.51	0.24	[0.10, 1.13]	NA	NA	NA	1.001	1088	1582
Additional Parameters
sigma	0.67	0.02	[0.64, 0.71]	NA	NA	NA	1.002	6374	2773

# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()

## Possible multicollinearity between b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.77). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

Device Based MVPA

Lognormal Model

formula <- bf(
  pa_obj ~ 
   
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_onlypushing", suffix))
)

## Warning: Rows containing NAs were excluded from the model.

pa_obj_log_digest <- digest::digest(pa_obj_log)

# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summarize_brms(
  pa_obj_log, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)

## Warning in summarize_brms(pa_obj_log, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.

	exp(Est.)	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	124.24***	7.52	[109.69, 139.58]	1.000	[0.94, 1.06]	0.000	1.002	967	1634
Within-Person Effects
Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing experienced	1.03	0.03	[ 0.97, 1.08]	0.807	[0.94, 1.06]	0.914	1.001	3221	2733
Daily pushing utilized (partner’s view)	1.02	0.02	[ 0.97, 1.06]	0.797	[0.94, 1.06]	0.979	1.000	5065	2748
Day	0.97	0.05	[ 0.87, 1.08]	0.718	[0.94, 1.06]	0.654	1.000	7350	3336
Actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily weartime	1.00**	0.00	[ 1.00, 1.00]	0.998	[0.94, 1.06]	1.000	1.001	4680	2898
Between-Person Effects
Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing experienced	0.97	0.06	[ 0.86, 1.10]	0.680	[0.94, 1.06]	0.635	1.001	2223	2645
Mean pushing utilized (partner’s view)	1.04	0.07	[ 0.92, 1.18]	0.739	[0.94, 1.06]	0.578	1.002	2058	2692
Mean weartime	1.00	0.00	[ 1.00, 1.00]	0.812	[0.94, 1.06]	1.000	1.000	3879	3246
Random Effects
sd(Intercept)	0.31	0.04	[0.24, 0.42]	NA	NA	NA	1.005	890	1790
sd(Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing experienced)	0.09	0.04	[0.02, 0.17]	NA	NA	NA	1.004	887	821
sd(Daily pushing utilized (partner’s view))	0.03	0.03	[0.00, 0.10]	NA	NA	NA	1.001	1398	1731
Additional Parameters
sigma	0.54	0.01	[0.52, 0.57]	NA	NA	NA	1.001	6040	3037

plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()

# Nothing significant, no plots

Affect

Gaussian

formula <- bf(
  aff ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan +
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss_onlypushing", suffix))
)

## Warning: Rows containing NAs were excluded from the model.

mood_gauss_digest <- digest::digest(mood_gauss)

summarize_brms(
  mood_gauss, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)

	Est.	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	3.78***	0.12	[ 3.56, 4.00]	1.000	[-0.11, 0.11]	0.000	1.013	560	1111
Within-Person Effects
Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing experienced	0.02	0.04	[-0.05, 0.08]	0.670	[-0.11, 0.11]	0.995	1.002	2823	2591
Daily pushing utilized (partner’s view)	0.09*	0.04	[ 0.01, 0.17]	0.984	[-0.11, 0.11]	0.715	1.002	1597	2207
Day	0.28***	0.09	[ 0.11, 0.46]	1.000	[-0.11, 0.11]	0.023	1.001	3871	3088
Actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing experienced	0.12	0.11	[-0.08, 0.32]	0.880	[-0.11, 0.11]	0.433	1.003	1269	1954
Mean pushing utilized (partner’s view)	0.08	0.11	[-0.12, 0.30]	0.786	[-0.11, 0.11]	0.554	1.002	1216	1958
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.61	0.08	[0.48, 0.79]	NA	NA	NA	1.011	648	1448
sd(Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing experienced)	0.07	0.04	[0.00, 0.16]	NA	NA	NA	1.000	1409	1126
sd(Daily pushing utilized (partner’s view))	0.14	0.04	[0.06, 0.24]	NA	NA	NA	1.002	1329	828
Additional Parameters
sigma	0.91	0.02	[0.87, 0.94]	NA	NA	NA	1.001	3887	2536

plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw

Reactance

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan +
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal_onlypushing", suffix))
)

## Warning: Rows containing NAs were excluded from the model.

reactance_ordinal_digest <- digest::digest(reactance_ordinal)

summarize_brms(
  reactance_ordinal, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)

	OR	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercepts
Intercept	NA	NA	NA	NA	NA	NA	NA	NA	NA
Intercept[1]	5.44***	1.57	[ 3.06, 10.06]	1.000	[0.84, 1.20]	0.000	1.001	3521	2882
Intercept[2]	9.89***	3.14	[ 5.45, 18.74]	1.000	[0.84, 1.20]	0.000	1.002	3598	2951
Intercept[3]	23.10***	7.89	[ 12.38, 46.79]	1.000	[0.84, 1.20]	0.000	1.001	3933	2899
Intercept[4]	81.88***	35.09	[ 38.53, 199.63]	1.000	[0.84, 1.20]	0.000	1.000	4569	3054
Intercept[5]	771.25***	625.85	[202.52, 4855.18]	1.000	[0.84, 1.20]	0.000	1.001	5246	2545
Within-Person Effects
Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing experienced	1.31*	0.14	[ 1.04, 1.62]	0.987	[0.84, 1.20]	0.214	1.001	3641	2884
Daily pushing utilized (partner’s view)	0.96	0.13	[ 0.75, 1.25]	0.615	[0.84, 1.20]	0.810	1.002	4044	2794
Day	1.60	0.73	[ 0.68, 3.95]	0.854	[0.84, 1.20]	0.182	1.002	4822	3245
Actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing experienced	2.32**	0.65	[ 1.35, 3.99]	1.000	[0.84, 1.20]	0.005	1.001	3394	2987
Mean pushing utilized (partner’s view)	1.25	0.35	[ 0.70, 2.16]	0.783	[0.84, 1.20]	0.365	1.001	3738	2562
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.85	0.21	[0.49, 1.34]	NA	NA	NA	1.001	1569	2204
sd(Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing experienced)	0.17	0.14	[0.01, 0.50]	NA	NA	NA	1.004	1241	1651
sd(Daily pushing utilized (partner’s view))	0.16	0.14	[0.01, 0.57]	NA	NA	NA	1.003	1626	1714
Additional Parameters
sigma	NA	NA	NA	NA	NA	NA	NA	NA	NA
disc	1.00	0.00	[1.00, 1.00]	NA	NA	NA	NA	NA	NA

plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()

## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.75), b_Intercept[4] and b_Intercept[3] (r = 0.82). This might lead
##   to inappropriate results. See 'Details' in '?rope'.

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan +
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance_onlypushing", suffix))
  #, file_refit = 'always'
)

## Warning: Rows containing NAs were excluded from the model.

is_reactance_digest <- digest::digest(is_reactance)

summarize_brms(
  is_reactance, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)

	OR	SE	95% CI	pd	ROPE	inside ROPE	Rhat	Bulk_ESS	Tail_ESS
Intercept	0.40**	0.11	[0.23, 0.68]	1.000	[0.83, 1.20]	0.002	1.000	3095	2830
Within-Person Effects
Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing experienced	1.37*	0.18	[1.06, 1.87]	0.993	[0.83, 1.20]	0.154	1.000	3238	2410
Daily pushing utilized (partner’s view)	1.16	0.26	[0.80, 2.12]	0.769	[0.83, 1.20]	0.509	1.000	2310	2170
Day	1.78	0.89	[0.69, 4.81]	0.874	[0.83, 1.20]	0.160	1.000	6984	2927
Actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Between-Person Effects
Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing experienced	5.39***	2.58	[2.38, 14.94]	1.000	[0.83, 1.20]	0.000	1.001	2293	2628
Mean pushing utilized (partner’s view)	2.37	1.15	[0.97, 6.82]	0.972	[0.83, 1.20]	0.057	1.001	2556	2640
Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	1.52	0.31	[1.03, 2.24]	NA	NA	NA	1.001	1431	2494
sd(Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing experienced)	0.18	0.15	[0.01, 0.56]	NA	NA	NA	1.003	1256	1876
sd(Daily pushing utilized (partner’s view))	0.42	0.27	[0.03, 1.15]	NA	NA	NA	1.002	1292	1879
Additional Parameters
sigma	NA	NA	NA	NA	NA	NA	NA	NA	NA

plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()

## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

conditional_spaghetti(
  is_reactance, 
  effects = c('pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_self_cw

Report All Models

process_model_component <- function(obj) {
  # Convert to string, modify, and evaluate
  name <- deparse(substitute(obj))
  if (report_hurdle) name <- paste0(name, '_hu')
  if (report_ordinal) name <- paste0(name, '_ordinal')
  return(get(name, envir = parent.frame()))
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  reactance_ordinal,
  is_reactance,
  
  stats_to_report = c('CI', 'pd'),
  
  model_rows_random = process_model_component(model_rows_random),
  model_rows_fixed = process_model_component(model_rows_fixed),
  model_rownames_random = process_model_component(model_rownames_random),
  model_rownames_fixed = process_model_component(model_rownames_fixed)
) 

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”

# pretty printing

summary_all_models <- all_models %>%
  print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),  
      "Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5), 
      "Mood Gaussian" = (length(all_models) / 5),
      "Reactance Ordinal" = (length(all_models) / 5),
      "Reactance Dichotome" = (length(all_models) / 5))
  )


export_xlsx(
  summary_all_models, 
  rows_to_pack = process_model_component(rows_to_pack),
  file.path("Output", "AllModelsOnlyPushing.xlsx"), 
  merge_option = 'header', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

  
print(summary_all_models)

	Subjective MVPA Hurdle Lognormal			Device-Based MVPA Log (Gaussian)			Mood Gaussian			Reactance Ordinal			Reactance Dichotome
	exp(Est.) pa_sub	95% CI pa_sub	pd pa_sub	exp(Est.) pa_obj_log	95% CI pa_obj_log	pd pa_obj_log	Est. mood_gauss	95% CI mood_gauss	pd mood_gauss	OR reactance_ordinal	95% CI reactance_ordinal	pd reactance_ordinal	OR is_reactance	95% CI is_reactance	pd is_reactance
Intercept	57.28***	[49.10, 66.97]	1.000	124.24***	[109.69, 139.58]	1.000	3.78***	[ 3.56, 4.00]	1.000	NA	NA	NA	0.40**	[0.23, 0.68]	1.000
Hurdle Intercept	3.79***	[ 2.27, 6.26]	1.000	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Conditional Within-Person Effects
Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily pushing experienced	0.97	[ 0.91, 1.03]	0.856	1.03	[ 0.97, 1.08]	0.807	0.02	[-0.05, 0.08]	0.670	1.31*	[ 1.04, 1.62]	0.987	1.37*	[1.06, 1.87]	0.993
Daily pushing utilized (partner’s view)	0.97	[ 0.90, 1.03]	0.847	1.02	[ 0.97, 1.06]	0.797	0.09*	[ 0.01, 0.17]	0.984	0.96	[ 0.75, 1.25]	0.615	1.16	[0.80, 2.12]	0.769
Day	0.91	[ 0.79, 1.06]	0.894	0.97	[ 0.87, 1.08]	0.718	0.28***	[ 0.11, 0.46]	1.000	1.60	[ 0.68, 3.95]	0.854	1.78	[0.69, 4.81]	0.874
Actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Daily weartime	NA	NA	NA	1.00**	[ 1.00, 1.00]	0.998	NA	NA	NA	NA	NA	NA	NA	NA	NA
Conditional Between-Person Effects
Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Mean pushing experienced	1.32*	[ 1.00, 1.71]	0.977	0.97	[ 0.86, 1.10]	0.680	0.12	[-0.08, 0.32]	0.880	2.32**	[ 1.35, 3.99]	1.000	5.39***	[2.38, 14.94]	1.000
Mean pushing utilized (partner’s view)	1.18	[ 0.90, 1.53]	0.883	1.04	[ 0.92, 1.18]	0.739	0.08	[-0.12, 0.30]	0.786	1.25	[ 0.70, 2.16]	0.783	2.37	[0.97, 6.82]	0.972
Mean weartime	NA	NA	NA	1.00	[ 1.00, 1.00]	0.812	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pushing experienced	1.55**	[ 1.17, 2.24]	0.998	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily pushing utilized (partner’s view)	1.84***	[ 1.38, 2.68]	1.000	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Day	0.61	[ 0.37, 1.00]	0.974	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Actionplan	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Daily weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean persuasion utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pressure experienced	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pressure utilized (partner’s view)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pushing experienced	0.46**	[ 0.26, 0.81]	0.997	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean pushing utilized (partner’s view)	0.62	[ 0.33, 1.08]	0.952	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Hu Mean weartime	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Random Effects
sd(Intercept)	0.37	[0.29, 0.50]	NA	0.31	[0.24, 0.42]	NA	0.61	[0.48, 0.79]	NA	0.85	[0.49, 1.34]	NA	1.52	[1.03, 2.24]	NA
sd(Hurdle Intercept)	1.25	[0.93, 1.74]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Daily pushing experienced)	0.06	[0.00, 0.16]	NA	0.09	[0.02, 0.17]	NA	0.07	[0.00, 0.16]	NA	0.17	[0.01, 0.50]	NA	0.18	[0.01, 0.56]	NA
sd(Daily pushing utilized (partner’s view))	0.09	[0.01, 0.18]	NA	0.03	[0.00, 0.10]	NA	0.14	[0.06, 0.24]	NA	0.16	[0.01, 0.57]	NA	0.42	[0.03, 1.15]	NA
sd(Hu Daily persuasion experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily persuasion utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pressure experienced)	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pressure utilized (partner’s view))	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pushing experienced)	0.46	[0.07, 1.02]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
sd(Hu Daily pushing utilized (partner’s view))	0.51	[0.10, 1.13]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
Additional Parameters
sigma	0.67	[0.64, 0.71]	NA	0.54	[0.52, 0.57]	NA	0.91	[0.87, 0.94]	NA	NA	NA	NA	NA	NA	NA

report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::report_packages()

• beepr (version 2.0; Bååth R, 2024)
• R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
• R.oo (version 1.27.0; Bengtsson H, 2003)
• R.utils (version 2.12.3; Bengtsson H, 2023)
• brms (version 2.22.0; Bürkner P, 2017)
• digest (version 0.6.37; Eddelbuettel D, 2024)
• Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
• bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
• lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
• DHARMa (version 0.4.7; Hartig F, 2024)
• wbCorr (version 0.1.22; Küng P, 2023)
• see (version 0.9.0; Lüdecke D et al., 2021)
• tibble (version 3.2.1; Müller K, Wickham H, 2023)
• R (version 4.4.1; R Core Team, 2024)
• openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
• ggplot2 (version 3.5.1; Wickham H, 2016)
• forcats (version 1.0.0; Wickham H, 2023)
• stringr (version 1.5.1; Wickham H, 2023)
• rvest (version 1.0.4; Wickham H, 2024)
• tidyverse (version 2.0.0; Wickham H et al., 2019)
• readxl (version 1.4.3; Wickham H, Bryan J, 2023)
• dplyr (version 1.1.4; Wickham H et al., 2023)
• purrr (version 1.0.2; Wickham H, Henry L, 2023)
• readr (version 2.1.5; Wickham H et al., 2024)
• xml2 (version 1.3.6; Wickham H et al., 2023)
• tidyr (version 1.3.1; Wickham H et al., 2024)
• knitr (version 1.49; Xie Y, 2024)
• kableExtra (version 1.4.0; Zhu H, 2024)

report::cite_packages()

• Bååth R (2024). beepr: Easily Play Notification Sounds on any Platform. R package version 2.0, https://github.com/rasmusab/beepr.
• Bengtsson H (2003). “The R.oo package - Object-Oriented Programming with References Using Standard R Code.” In Hornik K, Leisch F, Zeileis A (eds.), Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003). https://www.r-project.org/conferences/DSC-2003/Proceedings/Bengtsson.pdf.
• Bengtsson H (2003). “The R.oo package - Object-Oriented Programming with References Using Standard R Code.” In Hornik K, Leisch F, Zeileis A (eds.), Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003). https://www.r-project.org/conferences/DSC-2003/Proceedings/Bengtsson.pdf.
• Bengtsson H (2023). R.utils: Various Programming Utilities. R package version 2.12.3, https://github.com/HenrikBengtsson/R.utils, https://henrikbengtsson.github.io/R.utils/.
• Bürkner P (2017). “brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software, 80(1), 1-28. doi:10.18637/jss.v080.i01 https://doi.org/10.18637/jss.v080.i01. Bürkner P (2018). “Advanced Bayesian Multilevel Modeling with the R Package brms.” The R Journal, 10(1), 395-411. doi:10.32614/RJ-2018-017 https://doi.org/10.32614/RJ-2018-017. Bürkner P (2021). “Bayesian Item Response Modeling in R with brms and Stan.” Journal of Statistical Software, 100(5), 1-54. doi:10.18637/jss.v100.i05 https://doi.org/10.18637/jss.v100.i05.
• Eddelbuettel D (2024). digest: Create Compact Hash Digests of R Objects. R package version 0.6.37, https://dirk.eddelbuettel.com/code/digest.html, https://github.com/eddelbuettel/digest.
• Eddelbuettel D, Francois R, Allaire J, Ushey K, Kou Q, Russell N, Ucar I, Bates D, Chambers J (2024). Rcpp: Seamless R and C++ Integration. R package version 1.0.13-1, https://dirk.eddelbuettel.com/code/rcpp.html, https://github.com/RcppCore/Rcpp, https://www.rcpp.org. Eddelbuettel D, François R (2011). “Rcpp: Seamless R and C++ Integration.” Journal of Statistical Software, 40(8), 1-18. doi:10.18637/jss.v040.i08 https://doi.org/10.18637/jss.v040.i08. Eddelbuettel D (2013). Seamless R and C++ Integration with Rcpp. Springer, New York. doi:10.1007/978-1-4614-6868-4 https://doi.org/10.1007/978-1-4614-6868-4, ISBN 978-1-4614-6867-7. Eddelbuettel D, Balamuta J (2018). “Extending R with C++: A Brief Introduction to Rcpp.” The American Statistician, 72(1), 28-36. doi:10.1080/00031305.2017.1375990 https://doi.org/10.1080/00031305.2017.1375990.
• Gabry J, Mahr T (2024). “bayesplot: Plotting for Bayesian Models.” R package version 1.11.1, https://mc-stan.org/bayesplot/. Gabry J, Simpson D, Vehtari A, Betancourt M, Gelman A (2019). “Visualization in Bayesian workflow.” J. R. Stat. Soc. A, 182, 389-402. doi:10.1111/rssa.12378 https://doi.org/10.1111/rssa.12378.
• Grolemund G, Wickham H (2011). “Dates and Times Made Easy with lubridate.” Journal of Statistical Software, 40(3), 1-25. https://www.jstatsoft.org/v40/i03/.
• Hartig F (2024). DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models. R package version 0.4.7, http://florianhartig.github.io/DHARMa/.
• Küng P (2023). wbCorr: Bivariate Within- and Between-Cluster Correlations. University of Zürich. R package version 0.1.22. https://github.com/Pascal-Kueng/wbCorr.
• Lüdecke D, Patil I, Ben-Shachar M, Wiernik B, Waggoner P, Makowski D (2021). “see: An R Package for Visualizing Statistical Models.” Journal of Open Source Software, 6(64), 3393. doi:10.21105/joss.03393 https://doi.org/10.21105/joss.03393.
• Müller K, Wickham H (2023). tibble: Simple Data Frames. R package version 3.2.1, https://github.com/tidyverse/tibble, https://tibble.tidyverse.org/.
• R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.
• Schauberger P, Walker A (2024). openxlsx: Read, Write and Edit xlsx Files. R package version 4.2.7.1, https://github.com/ycphs/openxlsx, https://ycphs.github.io/openxlsx/index.html.
• Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, https://ggplot2.tidyverse.org.
• Wickham H (2023). forcats: Tools for Working with Categorical Variables (Factors). R package version 1.0.0, https://github.com/tidyverse/forcats, https://forcats.tidyverse.org/.
• Wickham H (2023). stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.5.1, https://github.com/tidyverse/stringr, https://stringr.tidyverse.org.
• Wickham H (2024). rvest: Easily Harvest (Scrape) Web Pages. R package version 1.0.4, https://github.com/tidyverse/rvest, https://rvest.tidyverse.org/.
• Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). “Welcome to the tidyverse.” Journal of Open Source Software, 4(43), 1686. doi:10.21105/joss.01686 https://doi.org/10.21105/joss.01686.
• Wickham H, Bryan J (2023). readxl: Read Excel Files. R package version 1.4.3, https://github.com/tidyverse/readxl, https://readxl.tidyverse.org.
• Wickham H, François R, Henry L, Müller K, Vaughan D (2023). dplyr: A Grammar of Data Manipulation. R package version 1.1.4, https://github.com/tidyverse/dplyr, https://dplyr.tidyverse.org.
• Wickham H, Henry L (2023). purrr: Functional Programming Tools. R package version 1.0.2, https://github.com/tidyverse/purrr, https://purrr.tidyverse.org/.
• Wickham H, Hester J, Bryan J (2024). readr: Read Rectangular Text Data. R package version 2.1.5, https://github.com/tidyverse/readr, https://readr.tidyverse.org.
• Wickham H, Hester J, Ooms J (2023). xml2: Parse XML. R package version 1.3.6, https://github.com/r-lib/xml2, https://xml2.r-lib.org/.
• Wickham H, Vaughan D, Girlich M (2024). tidyr: Tidy Messy Data. R package version 1.3.1, https://github.com/tidyverse/tidyr, https://tidyr.tidyverse.org.
• Xie Y (2024). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.49, https://yihui.org/knitr/. Xie Y (2015). Dynamic Documents with R and knitr, 2nd edition. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, https://yihui.org/knitr/. Xie Y (2014). “knitr: A Comprehensive Tool for Reproducible Research in R.” In Stodden V, Leisch F, Peng RD (eds.), Implementing Reproducible Computational Research. Chapman and Hall/CRC. ISBN 978-1466561595.
• Zhu H (2024). kableExtra: Construct Complex Table with ‘kable’ and Pipe Syntax. R package version 1.4.0, https://github.com/haozhu233/kableExtra, http://haozhu233.github.io/kableExtra/.